Successes and failures of Bethe Ansatz Density Functional Theory 
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The response of a one-dimensional fermion system is investigated using Density Functional Theory 
(DFT) within the Local Density Approximation (LDA), and compared with exact results. It is shown 
that DFT-LDA reproduces surprisingly well some of the characteristic features of the Luttinger 
liquid, namely the vanishing spectral weight of low energy particle-hole excitations, as well as the 
dispersion of the collective charge excitations. On the other hand, the approximation fails, even 
qualitatively, for quantities for which backscattering is important, i.e., those quantities which are 
crucial for an accurate description of transport. In particular, the Drude weight in the presence of 
a single impurity is discussed. 

PACS numbers: 71.10.Pm, 71. 15. Mb, 73.21. Hb 



I. INTRODUCTION 

Density Functional Theory (DFT) is the most efficient 
and powerful tool for determining the electronic struc- 
ture of solids. While originally developed for continuum 
electron systems with Coulomb interaction (is^ DFT has 
also been applied to lattice models, such as the Hubbard 
model.— i^i^i^ One goal of these studies was to develop new 
approaches to correlated electron systems: lattice mod- 
els often allow for exact solutions - either analytically 
or based on numerics - which hence can serve as bench- 
marks for assessing the quality of approximations that 
usually have to be made when using DFT. Very popular 
in solid state applications is the Local Density Approxi- 
mation (LDA) where the exchange-correlation energy of 
the inhomogeneous system under consideration is con- 
structed via a local approximation from the homogeneous 
electron system. Amongst others, LDA has been applied 
to study ultracold fermions in one-dimensional optical 
lattices^ Friedel oscillations in one-dimensional metals, * 
the Mott gap in the Hubbard model ^ and quantum spin 
chains 4 

For small systems and not too strong interactions, 
LDA in most cases produces reasonable results - which 
can be obtained with much less numerical effort than 
needed when using more accurate methods like exact 
diagonalization or the density matrix renormalization 
group. This led to the hope that the LDA may serve 
as a useful tool for large systems, where the numerical 
effort for the more accurate methods is too expensive. 

After recalling in the next section (Sect. II) the the- 
oretical background for the Bethe ansatz LDA, we will 
study in detail the LDA solution of spinless fermions in 
one dimension. In Sect. HI we determine the charge sus- 
ceptibility and discuss, in particular, questions of stabil- 
ity as well as the static and the dynamic response. Then 
(Sect. IV) we study the Drude weight (which can be re- 
lated to the conductivity) in the presence of a single im- 
purity, and present our conclusions in the final Sect. V. 
From the size-dependence of our results, wc find that for 
large systems, LDA predictions are qualitatively incor- 
rect even for weak interaction. 



II. FORMALISM 

We consider a one-dimensional model of spinless 
fermions described by the Hamiltonian 



•i+l 



'h.c.)-f (1) 



where cf{ci) creates (annihilates) a fermion at site i, t 
is the hopping parameter, V the nearest-neighbor inter- 
action, and Vi an arbitrary local potential. The lattice 
consists of L sites (the lattice constant a is set to one), 
and periodic boundary conditions are assumed. 

The lattice version'^ of DFT relies on the fact that there 
is a one-to-one correspondence between the potentials 
{vi} and the groundstate expectation values of the site 
occupations {n^}. Therefore it is - in principle - possible 
to express all quantities that can be obtained from the 
groundstate wave function as a function (or functional 
in the continuous case) of the densities. The site occu- 
pations as a function of the potentials can, of course, be 
found from derivatives of the groundstate energy with 
respect to the local potential. 
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On the other hand, in order to determine the potentials 
from the densities it is convenient to define the function 



F{{ni\) = min ( 



^!\T + V\^) 



(3) 



where ^E* — s- indicates that the minimization is con- 
strained to such wave functions 5* that yield the given 
site occupation, i.e., (\I'|ni|5') — rii. Here T and V are 
the kinetic and interaction part of the Hamiltonian ((T|) , 
respectively. The groundstate energy is obtained by min- 
imizing the function E{{ni\) = F{{ni\) -\- '^i^i with 
respect to n^. This yields the condition 



dF 







(4) 
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which, of course, is purely formal unless F or at least a 
reasonable approximation for it is available. 

A major step towards the practical implementation of 
DFT was the idea of Kohn and Sham^ to employ a non- 
interacting auxiliary Hamiltonian in order to calcu- 
late the groundstate density profile. In the present case 



(5) 



where the potentials vf have to be chosen such that in the 
groundstate of H'^ the site occupations rii are the same 
as in the interacting model. Performing the same steps 
as before, one obtains the conditions 

(6) 



On, 

Combining Eqs. ^ and 
d 



vt=0. 
] yields 



(7) 



where vf = V{ni+i + ni-i) is the Hartree potential, and 
wf^ is the so-called exchange-correlation potential. The 
minimization problem of DFT is thus mapped onto the 
diagonalization of the one-particle Hamiltonian sup- 
plemented with the self-consistency condition ([7]). How- 
ever, there remains the problem of finding a practical ap- 
proximation for the exchange-correlation potential vf^. 
Most DFT studies of lattice models have so far relied 
on the LDA where the groundstate energy density of 
the inhomogeneous system is approximated by the energy 
density of a homogeneous system at the same density. In 
the present case this quantity can be calculated from the 
Bethe ansatz equations; hence 

A 

dm 



J LDA 



(8) 



where (n) is the Bethe ansatz energy per site of a ho- 
mogeneous system with particle density n, and e^{n) the 
corresponding energy density in Hartree approximation. 



III. SUSCEPTIBILITY 

In order to assess the results based on LDA and to 
discuss their validity, it is appropriate to recall first the 
phase diagram of the model under consideration. In the 
homogeneous case and away from half filling one finds 
for all values of the interaction parameter V/t > — 2 a 
Luttinger liquid phase, i.e., there is no long range charge 
order and the low energy excitations are gapless collec- 
tive charge excitations. At half filling and for V/t > 2 the 
model exhibits long range charge order, and a charge gap 
opens. Figure 1 shows the exchange-correlation potential 
w™(n) obtained from Bethe ansatz, compare ([S]), for sev- 
eral values of the interaction strength V. Due to particle 
hole symmetry, we have v^^{l — n) = —v^^{n). Further- 
more, for V/t > 2 there is a discontinuity at n — 1/2, 
related to the opening of the charge gap. 
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FIG. 1: (Color online) Exchange correlation potential v^'^ of 
the one- dimensional spinless fermion model as function of the 
density n for several values of the nearest-neighbor interaction 
V. 
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FIG. 2: (Color online) Critical interaction strength Vc above 
which the LDA susceptibility is negative, indicating an insta- 
bility of the system with respect to charge ordering. In the 
infinite system, the stable region is localized near half filling, 
from Uc ~ 0.331 to 1 — Tie. For finite system size and weak 
interaction, LDA is stable for all densities. 



A. Stability 

First we study the stability of the homogeneous LDA 
solution by considering the charge susceptibility x(g). 
Generally, the susceptibility can be expressed as 



xiq) 



Xaiq) 



(9) 



where q is the wavevector, f^'^{q) the Fourier transform of 
f^^j — dvf^/drij, and V{q) — 2Fcosg; Xo{q) is the static 
susceptibility of the auxiliary system, given by {L —^ oo) 



xo{q) = 



1 



47rt sm{q/2) 



In 



sm(q/2) + sinkp 



sin(g/2) — sin fc/ 



(10) 



where kp is the Fermi wavevector. The stability bound- 
ary of the homogeneous density profile is determined by 
the condition that the static susceptibility becomes infi- 
nite and changes sign; this happens whenever the denom- 
inator in ^ vanishes, i.e., for V{q) + f^'^{q) = —Xo^iq)- 
Due to the logarithmic divergence of Xoiq) fo^' 9 ~^ 
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this is equivalent to the condition that V{2kp) + f^''{2kF) 
changes sign. Notice that within LDA the function 
f^'^io) ~^ /lda independent of q. Figure 2 shows the 
region of stabihty in the rt-y-planc obtained within LDA 
both for the infinite system, and for finite systems of 
length L = 100 and L = 1000, respectively. For L oo 
only systems with density near 1/2 and not too strong 
interaction are stable, further away from half filling the 
homogeneous solution is unstable for arbitrarily weak in- 
teraction. 

For an accurate determination of the phase boundary 
we investigate the weak interaction case in more detail. 
We find 



/] 



xc 

LDA 



BA 



eH) = -T/(2fcj.)+0(y2) (11) 



so that in first order in the interaction no conclusion 
about the stability can be drawn: the second order cor- 
rection to the groundstate energy, £2, is needed. Numeri- 
cally wc find that its second derivative with respect to the 
density, e'^iin), changes sign at ric ~ 0.331, thus limiting 
the range of stability to ric < n < 1 — ric a,t weak cou- 
pling. This result should be contrasted with the Hartree 
approximation [f^'^ = 0) where the homogeneous system 
is stable only for V{2kF) > 0, i.e., below quarter and 
above three quarter filling, and with the exact ground- 
state where a charge instability of the homogeneous sys- 
tem occurs only at half filling for V > 2t. 

There are, however, very pronounced finite-size ef- 
fects that strongly enlarge the actual region of stability 
within LDA. Since xoi'^kp) diverges only logarithmically 
with system size L, the critical interaction strength ap- 
proaches zero very slowly, Vc{L) ~ 1/Vln L. As a con- 
sequence, for finite systems and from weak to intermedi- 
ate interaction strength the homogeneous LDA solution 
is stable for all densities, as can be seen in Fig. 2 for 
L = 100 and L = 1000. 



B. Static response 




FIG. 3: (Color online) Static susceptibility xil) (in units of 
t"^) in LDA for V/t = 1 at half filling (L = 202, one-peak 
curve) and quarter filling (L = 204, two-peak curve). The 
symbols are results from exact diagonalization for systems of 
up to L = 20 sites. 





FIG. 4: (Color online) Imaginary part of the dynamical sus- 
ceptibility x^^°^{<l,i^) (in units of t"^). 



Here we investigate the static susceptibility, i.e., its q- 
dependence, in more detail. In Fig. 3 we show x^^^io) 
for V/t = 1 in comparison with the exact susceptibil- 
ity obtained from numerical diagonalization of small sys- 
tems. As to be expected, in the long wavelength limit, 
q ^ 0, perfect agreement is found. Technically, there 
is a cancellation between the susceptibility Xo^i^) ^ 
27r< sin fci? and the second derivative of the Hartree energy 
= —{2t/TT) sinkp + Vn? with respect to n = kp/n. 
Therefore 

which is the exact uniform susceptibility of the inter- 
acting system. Unfortunately, already the next to lead- 
ing contribution, ^ q^ , is not obtained correctly within 



LDA. At half filling the discrepancy between the LDA 
susceptibility and the exact one becomes more and more 
pronounced for q — > 2kp = tt. At q = 2kp the exact sus- 
ceptibility increases with the system size (not shown in 
the figure) and diverges with a power law, while in LDA 
there is only a cusp. The cusp value itself remains finite 
and approaches x'"'^^(7r) = 1.668/i for L ^ 00. 

At quarter filling x^^^(g) is very close to the exact 
susceptibility for q < 2kp, while for q > 2kp there is a 
clear discrepancy. For q — 2k p — 7r/2 the exact result 
again is strongly size-dependent and diverges for L — > 00, 
while within LDA the susceptibility diverges already at a 
finite systems size, since at quarter filling one is already 
outside the range of stability of LDA. 
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FIG. 5: (Color online) Charge velocity v as function of the 
interaction strength V for densities n = 0.1, 0.2, 0.3, 0.4, 0.5 
(from bottom to top). Exact values from Bethe ansatz (solid 
lines) in comparison with the results obtained within ALDA 
(symbols) from Eq. (114)1 . 
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FIG. 6: (Color online) Drude weight D for a half-filled system 
with a single impurity as function of the interaction strength 
V, for several values of the impurity strength Vimp/t = 0, 1, 
2, 3, 4, 5 (from top to bottom). The solid curve is the Bethe 
ansatz result for the clean system, the long- and short-dashed 
curves are LDA results for L — 102 and L = 202, respectively. 



C. Dynamic response 

DFT as presented in Sect. |TT] is a groundstate the- 
ory. However, a time-dependent generalization^^ of 
DFT is available which allows to calculate frequency- 
dependent response functionSfl? The dynamic suscepti- 
bility of the homogeneous system is of the same form 
as Eq. @ with the only differences that Xo{q) has to 
be replaced by xo(9, w), and f^^iq) by f^'^{q,uj). In 
a simple approximation, called adiabatic local density 
approximation^"^ (ALDA), one neglects this frequency de- 
pendence and uses the function f^'^ already known from 
LDA, f^^{q,Lj) — > /lda- The corresponding approxima- 
tion for the susceptibility is denoted x^^^^{q,i^)- 

Figure 4 shows the imaginary part of x^^^^iQ^^) iii 
the q-w-plane for a half- filled system and V/t = 1. A con- 
tinuum of excitations in the frequency range 2t \ sing| < 
uj < 4:tsm{q/2) is apparent, which can be identified 
with the particle-hole continuum. Note that the spec- 
tral weight of the particle-hole excitations vanishes for 
q ~* 0, a, feature - expected for a Luttinger liquid - which 
is reproduced surprisingly well in ALDA. 

Above this continuum we find a well-defined branch of 
collective excitations with linear dispersion for small q. 
Analytically, the dispersion of the collective excitations 
can be obtained from the singularities of the susceptibil- 
ity, i.e., by considering the zeros of the denominator of 
xiq, w). In the low frequency and small wavevector limit 
the susceptibility agrees with the Luttinger liquid result 



1 dN {qvf 
L dfi {qvY — uo'^ 



(13) 



where v is the velocity of the collective excitations. 
Within the adiabatic local density approximation the ve- 
locity is given by 



WALDA 




2V + n 



xc 

LDA 



■KVf 



(14) 



where vp — 2tsaikp is the Fermi velocity. To linear 
order in the interaction, walda agrees with the exact 
result. In Fig. 5 we compare walda with the exact 
value, obtained from Bethe ansatz, for various interac- 
tion strengths and densities. For small densities, there 
is nearly perfect agreement between the two values. The 
largest deviation occurs at half filling at the critical point 
V = 2t where the error is 2/^Jt: — 1, which is about 13 %. 
We want to emphasize that within the random phase ap- 
proximation, i.e., neglecting the ALDA correction factor 
f^'^ in Eq. p^ . one never obtains the correct charge ve- 
locity, except for V = Q. From the discrepancy between 
ualda a-nd the exact value one concludes that the fre- 
quency and wavevector dependent function f'^'^{q,uj) is 
non-analytic in the — > limit. 



IV. SINGLE IMPURITY 

Finally, we consider the spinless fermion model for a 
non-zero potential Vi. As an important example in this 
context, we consider the case of a single impurity, i.e., 
Vi = Wimp at the impurity site and Vi = elsewhere. It is 
well known that in one-dimensional interacting systems 
the reflection and transmission probabilities for scatter- 
ing at an impurity are strongly renormalizedi ^^i^^'^^ Here 
we address the question whether this renormalization is 
captured by the LDA - with negative conclusion. 

Let us first consider the simple picture for the origin of 
the renormalization that has been developed by Matveev 
et al.^'* Around an impurity, the density is disturbed, and 
in the presence of electron-electron interaction this mod- 
ulation (Friedel oscillations) leads to additional scatter- 
ing. In particular, the Friedel oscillations are the origin 
of enhanced backscattering. To linear order in the inter- 
action, the correction to the transmission probability for 
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a wavevector q close to kp is given bj^is 



ST = -2a%{l^%)\TL 



1 
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in agreement with the perturbative argument given in 
relation with Eq. 



(15) 



where To is the bare value; the dimensionless parame- 
ter a characterizes the interaction strength. It is given 
by the sum of a Hartree and an exchange contribu- 
tion, a = Q!H + ctxj with q;h — —V{2kp)/2'KVp and 
ctx = V{jd)l2'iTVp. By summation of the leading diver- 
gencies to all orders in the interaction using a rcnormal- 
ization group approach, it is found that even for a weak 
defect the transmission approaches zero as q ^ kp (re- 
pulsive interaction). Repeating the arguments leading 
to Eq. (jlSp within DFT and for a weak impurity, we 
find a -[V{2kp) + f'"'{2kp)]/2T:vp. Since /^Jj^ = 
—V{2kp), this singular correction to the transmission is 
zero, i.e., T is not renormalized in DFT-LDA. 

To substantiate this finding numerically we calculate 
the Drude weight for the single-impurity case. The Drude 
weight, D, is defined as the response of the system to a 
change of boundary conditions according to 



L d^E 
^2d^ 



(16) 



ip=0 



where E(ip) is the groundstate energy. The parameter 
if characterizes the twist in the boundary conditions: 
ip = corresponds to periodic, and ip = ir to antiperiodic 
boundary conditions fi^ii^ The Drude weight is closely re- 
lated to the transmission through the defect, and in the 
non-interacting system - where T is not renormalized - 
the size dependence of D is negligible. In the interact- 
ing system, on the other hand, the transmission coeffi- 
cient for (g — kp) « vp/L is relevant as discussed above. 
Correspondingly, the Drude weight increases (decreases) 
algebraically with system size for repulsive (attractive) 
interactionii^iiiii^ 

In Fig. 6 we present our LDA results for the Drude 
weight at half filling, for different system sizes {L = 102 
and L = 202) and different values of the impurity 
strength. Unlike the (numerically) exact results,-^ we do 
not observe any dependence on system size within LDA, 



V. SUMMARY 

We studied in detail the Bethe ansatz LDA for spinless 
fermions in one dimension. The numerical effort of the 
method is comparable to the Hartree (or Hartrec-Fock) 
approximation. A major improvement of LDA with re- 
spect to the Hartree approximation is the correct predic- 
tion of a non-charge-ordered groundstate for a large range 
of parameters. Both the static and the dynamic density 
response functions agree reasonably well with the exact 
results. In particular, for low density and q < 2kp, the 
static susceptibility obtained within LDA is almost indis- 
tinguishable from the exact one. In the dynamic case an 
impressive agreement for the velocity of collective charge 
excitations is found. 

On the other hand, the LDA does not capture the 
fact that the system is critical with respect to a charge- 
ordering phase transition. Whereas the exact suscepti- 
bility has a power law singularity at q = 2k p, the LDA 
susceptibility remains either finite or the system devel- 
ops charge ordering for very large system size. As a 
consequence physical quantities that are related to the 
2fcj7-periodic charge oscillations are described incorrectly 
within the local density approximation. Examples are 
the Friedel oscillations around a defect, the interaction- 
renormalization of the refiection and transmission prob- 
ability (and therefore the conductance) and the Drude 
weight. 

We conclude that for applications of DFT to 
one-dimensional systems improved exchange-correlation 
functionals are required. 
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